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Exact solution for low energy qnantum anharmonic vibrations in a long polymer chain 
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We propose the algorithm for determining quantum stationary states of periodic linear chain of 
atoms coupled by harmonic and third order anharmonic interactions (Fermi-Ulam-Pasta a problem) 
in the long wavelength limit within the resonant approach. These states can be encoded by sequences 
of integer numbers determining their energies and wavefunctions. Using these states we described 
the exact time evolution of a single phonon state showing coherent oscillations. The applications of 
theory to vibrational energy transport and quantum informatics are discussed. 
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Molecular vibrations determine the heat balance in 
nano-devices [l|, Q and can be manipulated similarly to 
electrons and photons and used to carry and process 
quantum information [^-Q. The molecules composed by 
chains of self-repeating monomers demonstrate an out¬ 
standing ability to transfer and convert energy because of 
delocalized normal modes (phonons) existing due to the 
translational invariance and propagating with the speed 
of sound as high as 10^ m/s in organic polymers [7Hlfll|. 

The increasing theoretical efforts have led to the re¬ 
markable progress in understandin g p olymer heat con¬ 
duction (see e. g. Reviews iEl lid . [l3 |) though the 
actual mechanisms of energy transport through anhar¬ 
monic oscillator chains remains unclear since its first nu¬ 
merical study by Fermi, Pasta and Ulam (FPU) [l^. The 
original numerical simulations of vibrational dynamics in 
the FPU a model including atoms coupled by harmonic 
and third order anharmonic interactions were targeted to 
reveal a thermalization accompanied by the loss of mem¬ 
ory about the initial state [15|. Instead a quasiperiodic 
dynamics has been found. This quasiperiodic dynam¬ 
ics has been interpreted introducing the solitary waves 
solutions and integrals of motion associated with them 
IS 14 1 3 12 1 suggesting that the system is integrable. 

It is natural to expect that the classical integrability 
should be reflected in the quantum mechanical proper¬ 
ties 18|. Quantum mechanical treatment is important 


for applications to real molecules since the thermal en¬ 
ergy ksT is usually smaller than the typ ical vibrational 
energy which is around 1000 cm“^ |10j (even at room 
temperature ksT ^ 200cm“^). Here we report the exact 
quantum mechanical solution for eigenstates and eigenen- 
ergies of the FPU a model in the lon g-wa velength limit 
restricted to the resonant interactions [l9l| . Below we in¬ 
troduce the quantum mechanical model, describe its so¬ 
lution and its application to the single phonon state time 
evolution. The details of the derivations of the solution 


are given in the Supplementary Materials [2l| . 

The normal modes (phonons) of the periodic (circu¬ 
lar) chain with the period a = 1 and length TN can be 
expressed as planar waves with the amplitude depending 
on the coordinate z = 1, ...N along the chain as x{z) = 


j^pN. In a periodic chain one has x{z A- N) = x{z) 
and the wavevector for each phonon can be expressed as 
q = 2TTn/N with an integer number n identifying each 
specific mode. In the long wavelength limit n N the 
phonon energy can be approximated by the linear dis¬ 
persion law En = hc\n\/N where c stands for the speed 
of sound (c = I in the FPU model mi). The system 
Hamiltonian in the harmonic approximation can be con¬ 
veniently expressed in terms of creation and annihilation 
operators b'l, bn for each mode n {Hq term in Eq. ([T])). 

The third order anharmonic interactions can be intro¬ 
duced using products of three 6-operators describing the 
phonon decay into two phonons blnbhbm+n or two phonon 
association backwards b\^j^j3rnbn conserving the total 
wavevector within the long-wavelength (low energy) limit 
(cf. [3 mi) 22|). We restrict the consideration to the 
only fully resonant interactions expressed by the terms 
with both parameters m and n either positive or nega¬ 
tive [3|. Indeed, the harmonic energy does not change 
in such process {h{\m -I- n| — \m\ — |n|)/A^ = 0) while it 
changes by a large harmonic energy 26min(|TO|, |n|)/A^ 
in the opposite case. Within this resonant approach the 
system Hamiltonian can be separated into two parts as¬ 
sociated with positive and negative wavevectors. The 
positive wavevector part can be expressed as [3, El (see 
also Ref. 2l|, Sec. I; we set the Planck constant 6, = 1, 
the negative wavevector part can be studied similarly) 


Hres — 1^0 


y/2N^ 


V, ^0 = + 1 / 2 ), 


n>0 


E = 2 X Vmnim + n) (6+6+6„+„ -b 6++„6„6„) .(1) 

m,n>0 

Here the parameter a is the relative anharmonic interac¬ 
tion used in Ref. [3|. 

Each eigenstate of Eq. o is determined by a super po¬ 
sition of multiphonon states defined by population num¬ 
ber sequences {v} = pi, 1 / 2 , ■■■Vn) Pi = ^l^i) as 
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where c{{iy}) are modified wavefunction amplitudes for 
each state {i/}. After this modification adding factors 
iV ({:/}) the action of the Hamiltonian is expressed by 
integer numbers Eq. 

Since the harmonic energy is conserved in the reso¬ 
nant approximation {[Hres,Ho] = 0, see Eq. (P)) the 
problem is reduced to the diagonalization of the dimen¬ 
sionless anharmonic interaction Hamiltonian V and each 
contributing multiphonon state should possess the same 
harmonic energy n/N satisfying the identity 

n 

= n. (3) 


Then each state can be characterized by its principal 
quantum number n. Each selection of population num¬ 
bers {vi} satisfying Eq. ([3|) corresponds to a certain inte¬ 
ger partition of the number n [23| representing a way of 
writing n as a sum of positive integers. For instance for 
the principal quantum number n = 3 there exist three 
different partitions (3 = 1-|-1-|-1, 3 = 2-|-l, 3 = 3, see 
Fig. P corresponding to the phonon population numbers 
= (3,0,0), (1,1,0) and (0,0,1). 

Here we propose the algorithm to determine eigen¬ 
states of the dimensionless Hamiltonian V Eq. o and 
corresponding eigenenergies. One can describe the possi¬ 
ble eigenstate using the sequence of p-1-1 integer numbers 
{k} = (ko,ki, ...kp), such that ko = n and kp = 0 {n is 
the principal quantum number). The following rules de¬ 
termine the eigenstate and energy corresponding to this 
sequence. 

1. Short sequences (n,0) correspond to the stationary 

states will all amplitudes equal to one {c{{iy}) = 1, see 
Eq. @) and energies e(n,o) = ~ l)/2- This can be 

proved using direct substitution (cf. Eq. (ITOH . 

2. The eigenstate ({;/}) (if non-trivial) correspond¬ 
ing to the given sequence {fc} = {ko,ki, ...kp) can be 
defined using the eigenstate C{fc_}({m}) for the reduced 
sequence {k-} = {ki,...kp) with the principal quantum 
number ki (^^ anria = ki) obtained removing the zeroth 
term from the original sequence. The connection between 
two solutions and their energies can be expressed as 


{m} 

n{n — 1) ki {ki — 1) 

Cf/c} = e{fe-} H- 2 --2-' 


( 4 ) 


where the summation is taken over all partitions {m} of 
the number ki and the functions '0{m}({i^}) are given by 
the products of associated Laguerre polynomials 


i=i 


( 5 ) 


To obtain the eigenstate described by the sequence {k} 
this algorithm should be repeated p times beginning with 


the sequence (fcp_i,0) corresponding to the all ones so¬ 
lution. 

3. If the sequences {k} are chosen strictly decreas¬ 
ing and satisfying the rule ki-i — ki > ki — ki+i then the 
number of sequences is equal to the number of partitions. 
The numerical studies show that up to the maximum ac¬ 
cessed principal quantum number n = 25 the eigenstates 
generated using this sequences following the above algo¬ 
rithm form the complete basis of eigenstates all normal¬ 
ized by one and orthogonal to each other. Unfortunately, 
we cannot give a general proof of this statement for ar¬ 
bitrarily n though the normalization by 1 is proved for 
some groups of generated states [21| Secs. H, V. The use 
of basis functions with smaller quantum numbers ki < n 
to describe the partitions of larger number n does not 
conflict with the completeness of the basis because the 
populations numbers are dependent of each other, i. e. 
they are bound by the “harmonic energy conservation 
law” Eq. @. 

The dimensionless energy of the eigenstate obtained 
repeating p — 1 times the iteration procedure Eq. (jj]) for 
a certain sequence {fc} is given by 
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FIG. 1: Partitions, eigenstates and dimensionless anharmonic 
eigenenergies for principal quantum numbers n = 1,2, and 3 
(See Ref. [2j|, Sec. VI for detail). 


Using these eigenstates one can describe the exact time 
evolution of the single phonon state assuming that at 
time t = 0 there was only one phonon in the harmonic 
state with the energy n/N (cf. Eq. ([T])). It can be shown 
(see below Eq. (fT^ f that the probability that the system 
remains in the single phonon state oscillates with the time 
as (see Fig. HD 


Pu{t) = 


sin^ 
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( 7 ) 


These oscillations and the oscillation period dependence 
on the anharmonic interaction and the system size are 
similar to the behaviors discovered in Ref. [l^ . 

To illustrate the proposed algorithm we show in Fig. [1] 
all eigenstates for n = 1, 2 and 3 obtained using strictly 
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decreasin g se quences of quantum numbers k (for detail 
see Ref. [2l|, Sec. VI). Basis states divided by corre¬ 
sponding normalization factors (see Eq. @) are repre- 
sented by diagrams [ 2 ^. For n = 1, 2 the strictly decreas¬ 
ing sequences represent the complete basis of eigenstates 
of the problem. For n = 3 three sequences (3,0), (3,1, 0) 
and (3, 2,1,0) represent the complete basis set in agree¬ 
ment with the proposed algorithm. The sequence (3,2,0) 
leads to the zero wavefunction and there are many more 
of such sequences for larger n among strictly decreasing 
sequences. This is because the total number 2" of strictly 
decreasing sequences {k} is much greater than the total 
number of partitions depending on n as [ 2 ^ 


p{n) 
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FIG. 2: Time evolution of the probability for the system to 
remain in its initial single phonon state with the principal 
quantum number n (see Eq. 0). 

Consider the derivation of the results, described above. 
To find eigenstates of Eq. © we can use the wavefunc¬ 
tion in the form Eq. @ for some principal quantum 
number n expressing the harmonic part of energy. The 
Schrodinger equation for the dimensionless anharmonic 
Hamiltonian V (cf. Eq. ([IJ) can be written as 

a,b 

+ + (9) 

a.b 

The introduced a?-operators raise or lower the popu¬ 
lation index Va of the amplitude by 1 (for example 
for n = 3 one can express the wavefunction amplitudes 


as ,i.' 2 , 1^3 and then X 2 

Population numbers cannot be negative; fortunately the 
related terms disappear in Eq. because of the zero 
factor I'ab'b so there is no need to care about them. 

We begin with the consideration of the solution de¬ 
termined by the sequence (n, 0 ) which is claimed to be 
C{^} = 1 for all partitions {iz} and it should have the 
eigenenergy e = n(n — I)/2 Eq. ([5]). Indeed, assuming 
all identical amplitudes one can rewrite Eq. @ for 
some specific partition {i/} as 
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1 
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abi^ai^b 

a,b 



( 10 ) 


Since for each partition {v} in Eq. ([2]) one has ava = 
n (see Eq. ([3])) we got e = n{n — l)/2 as in Eq. ®. 

The recursive algorithm of finding solutions described 
above is based on the following property of solutions of 
Eq. ([9]). If the set of amplitudes dy^} for the principal 
quantum number m describes the solution with the en¬ 
ergy Cd it can be used to generate another solution cy^y 
(if it is non-trivial) with the principal quantum number 
n and eigenenergy Cc defined using Eq. ® as 


fm} 

n(n — 1) mim — 1) 

Cc = Cd H-- nm H---. 


( 11 ) 


This result can be derived substituting the solution in the 
form of Eq. (HU into the Schrodinger equation Eq. ® 
and simplifying it using the properties of Laguerre poly¬ 
nomials [ 2 ^ (see Sec. Ill in Ref. 211 for detail). The 
algorithm Eq. m of finding eigenstates and eigenener- 
gies is implemented in Matlab codes which are the part 
of Ref. [21| (Sec. VII in the text and the code files) and 
the reader is strongly encouraged to use them and Sec. 
VI in Ref. [2l| to verify the relevance of the proposed 
algorithm. 

Choosing kg = n, ki = m one can identify the expres¬ 
sion of amplitudes cy^y in terms of amplitudes dy^y as 
the first iteration step Eq. O which can be repeated 
(backwards) arbitrarily number of times until the termi¬ 
nation at kp = 0 where all amplitudes should be set equal 
unity as described above. The corresponding evolution of 
energy by the set of shifts Eq. (HU leads to Eq. (|6]). Thus 
each integer number sequence defines the eigenstate and 
eigenenergy of the problem if this eigenstate is nontrivial. 

The amplitudes cy^yHv}) for the specific system eigen¬ 
state determined by the sequence {fc} of p -I- 1 quantum 
numbers and taken for the specific population number set 
{i/} can be expressed by means of the generating func¬ 
tion calculated using Eq. (jd]) in the form (see Ref. [21|, 
Sec. IV) 
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G{{v},{y}) = ^C{k}{{v})kfkf 

{fe} f,gii<f<g<p) 




n [l + 2/f + -{yi-yp-iy 

i=i 


( 12 ) 


Wavefunction amplitudes for the sequence {k} are ex¬ 
pressed by the coefficients of the generating func¬ 
tion expansion with the term containing the product 

fci k'2 — 1 

Vi y2—yp-i ■ 

Eq. (fT^ can be used to define all states having non¬ 
zero overlap with the single phonon state \n > repre¬ 
sented by the population number sequence Vk = ^kn, 
which is needed to describe the time evolution of this 
state Eq. ©• We consider only strictly decreasing se¬ 
quences {A:} which is sufficient to get a complete expan¬ 
sion of the initial state over eigenstates of the problem 
as shown below. Then for the state Vk = Skn one can 
leave only unity from the second factor in the right hand 
side of Eq. m because otherwise the power of some 
of variables yi will approach or exceed n which is not 
acceptable since these powers should form strictly de¬ 
creasing sequence n > ki > k 2 > k^... > kp = 0. 
The only acceptable choice of the contribution from 
the first product satisfying the sequence selection re¬ 
quirements can be made taking the composite product 
yi X ( 2 / 12 / 2 ) X ( 2 / 12 / 22 / 3 )--- X ( 2 /i 2 / 2 -- 2 /p-i)- The related se¬ 
quence is given by {n,p — l,p—2,p— 3, ...0) and it deter¬ 
mines the stationary state 'kp characterized by the energy 
ep = {—n{n — l)/2 + n{n—p)) (cf. Eq. ([6|)). There are n 
such sequences and associated eigenstates 'kp determined 
by the integer number p changing from 1 to n (cf. the 
solutions for eigenstates in Ref. [2l|, Sec. VI for n < 5). 
In all cases (see Ref. [2l|, Sec. V) the absolute value of 
the amplitude of the wavefunction in the single phonon 
state is equal unity, and it can be recalculated as l/\/u 
for the normalized by 1 representation Eq. ©■ Thus we 
found n contributing eigenstates having the overlap inte¬ 
gral Cp = ^l\pn with the single phonon state of interest. 
The expansion of the single phonon state over the basis 
of those n states is complete since nCp = 1. 

Assume that at time t = 0 the system is in a single 
phonon state \n >. Then the time evolution of the prob¬ 
ability Pn{t) to find the system in this specific state can 
be evaluated expanding this state over the previously es¬ 
tablished set of n eigenstates 'kp as 


Pn{t) = 


P =1 


< 4'p|n > 




Since all overlap matrix elements are identical so that 
I < 'kp|n > P = 1/n the sum in Eq. (fT^ takes the form 
of the geometric series X)p=i The evaluation of 

this geometric series results in Eq. ©. 


Consider some properties of eigenstates of the problem. 
Since Eq. 0 has a symmetry with respect to the trans¬ 
formation e —>• —e, Cjj,} —>■ Cjj,} • (—all eigenstates 
with non-zero energy enter in pairs (this conclusion is il¬ 
lustrated in Ref. 2l| for n < 5). Particularly the all ones 
state determined by the sequence (n, 0 ) possessing the en¬ 
ergy n(n— l )/2 has a corresponding state determined by 
the sequence (n,n — l,u — 2 ... 0 ) possessing the opposite 
energy (see Eq. d®]), Ref. [HI). The first states pos¬ 
sesses the maximum of energy because it has all positive 
amplitudes, while the anharmonic Hamiltonian has all 
positive matrix elements. Consequently the second state 
possesses the energy minimum and energies of all other 
states belong to the domain (—n(n — l)/ 2 ,n(n — l)/ 2 ). 
Since all energies are expected to be expressed by integer 
numbers Eq. and the number of states (partitions, 
Eq. ([ 8 ])) grows with the principal quantum number n 
faster than any power of n the strong degeneracy is ex¬ 
pected at large n, reflecting the integrability of the prob¬ 
lem. 

The fourth order anharmonic interaction can be in¬ 
troduced within the resonant approach similarly to Eq. 
©• It will represent the /3 FPU model characterized by 
some constant /3 of the fourth order anharmonic interac¬ 
tion. The preliminary numerical study of this problem 
does not lead to the analytical solution; yet, the small 
modification of the resonant /3 model by adding to the 
original Hamiltonian the diagonal term in phonon popu¬ 
lation numbers proportional to the expression 


E 


a Va{l + Va) - 3 “ 


12 


E 


a^Vn 


(14) 


makes the problem eigenstates identical to those of the 
resonant FPU a model. The accurate analysis of the /3 
FPU problem will be performed separately. 

In addition to vibrational energy transport the ob¬ 
tained solution can be relevant for quantum informat¬ 
ics because of its connection to the number theory (^ . 
Therefore its realizations employing interacting Joseph- 

251 is of interest and the 


son junctions or cold atoms 
model Eq. © can be hopefully implemented there with 
a high accuracy. 

In realistic polymers the breakdown of integrable be¬ 
havior and transition to chaos are possible due to the 
omitted “non-trivial resonances” (due to high order pro¬ 
cesses messing up positive and negative wavevectors) in¬ 
evitably leading to the ergodic behavior in a classical 
system according to Ref. [l9|. We hope that in a quan- 
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turn system an integrable state can be stable because of 
the phase space discreteness jl^. This problem needs a 
separate study. 

Even at small phonon energy n ^ 1 the deviation of 
the phonon dispersion law from the linear one 1 /N^ Ref. 
[2lj, Sec. I can destroy Fermi resonances when it exceeds 
the Fermi resonance amplitude a/N"^ Eq. Ilj. Conse¬ 
quently, the proposed theory can be applicable for a suf¬ 
ficiently large system size TV > 1/a. For organic polymers 
a ^ 0.1, and the regime N > 1/a ^ 10 is quite accessi¬ 
ble. The long wavelength limit requires a typical phonon 
energy jN {n is a principal quantum number and the 
typical phonon energy is taken as a thermal energy cor¬ 
responding to the total energy n/N [i^) to be less than 
the Debye energy, which is of order of 1 within the FPU 
model. More accurate analysis of applicability limits for 
the present solution requires separate investigation. The 
obtained solution should not be very sensitive to defects 
in the long-wavelength low energy limit 11, 12, 27, 2^ . 

Thus the family of analytical solutions is found for 
eigenstates of the quantum mechanical many-body prob¬ 
lem of one dimensional acoustic vibrations coupled by 
the resonant anharmonic interactions. According to the 
numerical studies this family forms the complete set of 
eigenstates. Eigenstates are described by sequences of 
quantum numbers possibly representing the integrals of 
motion of unclear nature which calls for further theoreti¬ 
cal studies. Practically the present model on a one hand 
is closely related to the vibrational energy transport in 
molecular chains and on the other hand it is connected 
to the number theory thus having a potential interest 
in quantum informatics. Therefore we hope that this 
work will stimulate experimental efforts to implement 
the present model using cold atoms and/or Josephson 
junctions and contribute to understanding the thermal 
conductivity of polymers. 

This work is supported by the National Science Foun¬ 
dation (CHE-1462075). 
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Supplementary Materials 


DERIVATION OF THE RESONANT ANHARMONIC HAMILTONIAN 


We begin with the derivation of the resonant Hamiltonian for the FPU a problem describing anharmonic vibrations 
of atoms forming periodic chain. Remember that a problem involves the third order anharmonic interactions only. 
The system Hamiltonian is defined following the seminal work of Fermi, Pasta and Ulam 15| as 


^ 1 ^ 1 


2=1 

N-1 


+ w + w(a;Ar - xi) 


2=1 


(15) 


This Hamiltonian includes harmonic (first term) and third order anharmonic (second term) interactions. 

Normal modes diagonalizing the harmonic part of the Hamiltonian Eq. (USD can be introduced as periodic waves 
or phonons (assuming even N without the lack of generality) 




-iqn 

c ‘^n 5 


q n—1 

27rm N N 

Similar transformation should be performed for atomic momenta. 

Using these new coordinates and momenta one can reexpress the Hamiltonian Eq. (USD as 

1 


(16) 


-H = ^ E {P<lP-<l + ^luqU-q) + 


9>0 


+ E ^ 91 + 92 + 93 ( 1 -e*'^')(l-e*®)(l-e”®)u 9 iU,,M, 3 , 


uj{q) = 2 I sin(g/2) |, 


A 


91+92+93 


= E 


+ 7722 + 2723,nAT 5 


(17) 


where ui{q) represents the vibrational frequency of a phonon with the wavevector g, the 5 symbol, 5m,n = 1 for m = n 
and 0 otherwise, stands for the Kronecker symbol and represents the integer quantization number corresponding 
to the wavevector Eq. (USD- Since at g <C 1 one has w ~ g the speed of sound is equal to unity within the FPU a 
model. 

It is convenient to express the Hamiltonian in terms of creation and annihilation operators of vibrational modes 
defined as 


Un = 


y^( 6 + + &_,), p, = - b.,). 


Consequently one can represent the harmonic Hamiltonian in its standard diagonal form 

Ho = ^hhjg{b+bq + 1/2), 


(18) 


(19) 


while anharmonic interactions can be expressed as 


Is = 1 ^91+92+936 

^ ^ 91,92.93 

ah^ 


+i 


3+V 


E A 

91,92,93 


.- 91+<12+<13 sin(gi/2) sin(g2/2) sin(g3/2) 

+1 sin(gi/2) sin(g 2 / 2 ) sin(g 3 / 2 )| '' 

.- 91+92+93 sin(gi/2) sin(g 2 / 2 ) sin(g 3 / 2 ) 


91+92+93® 


+1 sin(gi/2) sin(g 2 / 2 ) sin(g 3 / 2 )| 


{bt,b+^-q3-btg^bg,bg,) + 

{K1K2K3 - Kbqibq^) ■ 


( 20 ) 
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Assuming that the only low energy vibrations are considered gi, 92-93 ^ 1 (i- e. ni, 712 - 71,3 <C -/V ) we can restrict the 
summation in the A symbol definition in Eq. (II3 to the case n = 0 assuming the conservation of quasi-momentum 
and set sin{q/2) k, q/2 everywhere in Eqs. (fTOl) and (1^01) . Then the Hamiltonian takes the approximate form 


H—— ^ I 71 I (6 +-I-1/2)-I- 


a/i 2 


—1 


mn{m + n) / + + + . 

/ , n - 7 - - - ^ \^m^n^m+n ^m+nPrnbn) 


^ \/\mn{m + n)\ 


ah^ 


^ mn(m + n) . + + + . 


3\/8A2 ^ y'|7T17l(771-k7l)| 


( 21 ) 


Finally we leave only “fully” resonant processes conserving both energy and quasimomenta and perform the unitary 
transformation 6^ —>■ , bn —>■ —ibn- Then the Hamiltonian can be separated into positive and negative wavevector 
71,771 parts which can be treated separately. The positive wavevector part can be written as 


Hres =Ho + Ho = -^Y. ^(^nbn + 1/2); 

^ n>0 

^ = 'Jmn{m + n) (5+6+&™+„ -h b1;^j^J)mbn) - 

771, n>0 


( 22 ) 


This Hamiltonian is studied within the main body of the manuscript. Since the harmonic energy is conserved (the 
harmonic Hamiltonian Hq commutes with the whole Hamiltonian) the problem can be reduced to the diagonalization 
of the anharmonic dimensionless Hamiltonian V. 

The correction to the linear dispersion law for the phonon frequency Eq. 03 for the given quantum number n can 
be estimated as Suj ^ (n/N)^. This correction is always smaller than the third order anharmonic interaction Eq. (1211) 
at sufficiently large system size N because the latter interaction decreases with this size as N~‘^. 


NORMALIZATION OF THE SOLUTION FOR THE SEQUENCES (n, 0). 

For each partition {v} of the number 71 representing the corresponding sequence of population numbers the squared 
amplitude of the wavefunction for the state determined by the sequence (71, 0 ) is given by (we assumed all 
amplitudes = 1) 






1 


n- 




(23) 


One can calculate the sum of squared amplitudes over all possible partitions using the generating function G{x) 
defined as 


= (24) 

where the summation is over all possible partitions {v} for all integer numbers 71. Then this sum can be represented 
by the power series G{x) = where the coefficients Nk are defined as the normalization factors for each 

principal quantum number k, i. e. 


Gyk 

Using the definition of the amplitudes Eq. (ESI) one can evaluate the generating function as 


G{x) = \ay^y ^xp ( T " 


lO 


1 — X 


(25) 


(26) 
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Since all coefficients in the expansion of Eq. (|26)) with respect to the powers of x are equal to 1 we can conclude that 
all eigenstates (n, 0) are normalized by 1 for arbitrarily n. 

Using a similar method one can also prove the normalization by one for the quantum number sequence (n, m,m — 
1, ...0) with m < n (see Sec. Numerical probes show that the normalization by one holds for all non-trivial 
eigenstates generated from sequences determined by the algorithm formulated within the main text. Yet we cannot 
prove this statement in a general case. 


DERIVATION OF THE RECURSIVE EQUATION 


The original Schrodinger equation for the modified wavefunction amplitudes 
respect to the anharmonic Hamiltonian V has the form 

a.b 

a,b 


(Eq- (HSI) with 


(27) 


where the operator h describes the action of the anharmonic Hamiltonian V in the representation of modified am¬ 
plitudes Remember that ^-operators x^ raise or lower the population index Va of the amplitude c{{i'}) by 

1 . 

Assume that the amplitudes are expanded in terms of the basis set composed by Laguerre polynomial products 




i=i 


(28) 


as 


C{u} (29) 

{m} 

where the Kronecker symbol Sj 2 .wi,n is equal to one for the population number set {vi} satisfying the conservation 
law 

n 

ivi = n. (30) 

i=l 

This symbol defines the all unities wave function corresponding to the sequence (n,0). 

We are going to show that the amplitudes dym] in Eq. (Hi) can be chosen in the way that they differ from zero 
only for sequences {to} representing integer partitions of some number m suggesting 

irrii = TO. (31) 


Moreover these amplitudes satisfy the equation 


1 

2 


y^^abvgivi, 

a.b 


n(n — 1) m(m — 1) 

e---h nm --- 


^{m} — 


dab)ya Vb yl+bd{m} + 


a.b 


(32) 


where raising and lowering operators act on the indices to. This equation is almost identical to Eq. (HZl) except 
for the redefinition of energy that determines the recursive procedure defining the energy of each specific eigenstate 
in terms of the set of the associated quantum numbers. 

To derive this equation one can seek the solution in the form 




( 33 ) 
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where the sum is taken over all integer partitions of a certain number m Eq. (ED). 

The wavefunction amplitudes Eq. (l2^ can be expressed as the results of the action of operators iIj = V^{m}({«}) on 
the unit wavefunction ipj Eq. ([29]). One can then represent the action of the Hamiltonian in Eq. ED as 

+ [h, (34) 

It was shown in the main body of the manuscript that the state ^/>/ with all amplitudes c equal 1 is the eigenstate of 
the problem with the eigenenergy n{n — l)/2 so that hipj = n{n — l)/2^/;/. Consider the commutator [h,ijj] related 
part of the problem. 

Since we are interested in the Laguerre polynomial dependence of population numbers n it is convenient to introduce 
the different notations for them as 




(35) 


which is easier to follow during the derivation below. 

The commutation rules that can be used to evaluate the expressions in Eq. (1341) can be summarized as following 


[Xa , P^A’^c)] = 6ac{P^S^c “ 1 ) - P^S^c))Xa = -SacP^^-l{l'c “ l)x^ ; 

[^6 > [Xa 1 PmS^c)]] — 5ac^bcPm^-2{^c — 2 ); 

[xt.P?n,(yc)] = 6ac(.PmS^c + 1 ) - P^^{Vc))xt = 5acP^n,-l(Vc)xt ] 

[xt, [xt, P^^i'^c)]] = 6acSbcP^,-2{’^c)- (36) 


These rules are based on identities for Laguerre polynomials 2^ that have been used during the derivation of com¬ 
mutators in Eq. (l36l) . These identities can be rewritten in the notations of Eq. (|35p in the form 


p^uM) - p^uM - 1 ) = - 1 ), 

P^.-M - 1 ) - - 2 ) = - 2 ); 

P^^-M + 1) - P“_i(^c) = P^_2(^^c). (37) 

Eq. (1361) permits us to bring all raising or lowering x operators to the right hand side to act directly on the 
Kronecker symbol wavefunction tjjj = ivi,n as 


ahvai^b - Sab)Xa a?;, x'^+bi’i = abua{vb - Sab)'ipl 


(38) 


or 


(a -I- b)va+bXaX'^x^_^_^4’i = {a + b)i^a+b4’i- (39) 

The action of the product of three operators does not change the wavefunction if 7 ^ 0 in Eq. (|551) or if 7 ^ 0 
in Eq. (IMl) bringing it to zero otherwise. However in that case (e. g. Va = 0 or r';, = 0 in Eq. (1551) 1 the operator 
action results in the zero answer because of the presence of operators Va, ^b or Va+b- Consequently we can skip the 
product of x-operators if they are applied directly to the unit function tpi and then we can also skip the unit function 
■ 0 / to avoid the complexity in the notations. 

The main target of the further mathematical consideration is to reexpress all polynomials bringing their shifted 
arguments Vi ± 1 back to Vi by means of shifting the parameters rrii and then sim plify ing the like terms that should 
lead to the target Eq. (l32)) . Additional identities related to Laguerre polynomials [2^ will also be used used for this 
purpose. These identities are summarized below 


’^cP^,-i{i^c - 1 ) = mcP^^{vc) + P^^_i(j^c)/c; 

l/e(j^c- l)P^„-2(j^c- 2) = mc(TOc- l)P^^(j^c) + 2 (TOc- l)P^,-l(j^c)/c-kP^^_2(j^c)/c^ 

Cl'cP^.-lil'c) = CmcP^^{Vc) + {cmc + 1 - c)P^^_j^{Vc) + P^„_2(j^c); 
ci'cP^S'^c) = c{mc + l)P^^+i{i'c) + {crUc + l)P^^{i^c) + P^,-i{i'c)- (40) 

Also the identity defining the principal quantum number n Eq. (130]) will be used. 

To evaluate the commutator of the Hamiltonian with the arbitrarily product of Laguerre polynomials 0{m}({^}) 
we consider the most general form of the commutators with the first (top line) and second (bottom line) parts of 
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the Hamiltonian Eq. (EH) enumerated by the numbers (1) for the first part with a b, (2) for the first part with 
a = b, (3) for the second part with a ^ b, (4) for the second part with a = b. Without restricting the generality of 
the consideration we assume a > 6 in the cases a ^ b. The Laguerre polynomials commuting with the part of the 
Hamiltonian under consideration are skipped for the sake of simplicity 

One can evaluate the commutators breaking it into more enumerated parts until the desirable form of expressions 
designated by the symbol is attained for each part. For the first type of commutators enumerated by the index 
(1) we get 


(1) = avahvb[x^ = 

= Jz., - l)P^,(z/6 - l)[x++„P“t',(z^a+6)](1.2). (41) 


The first term (1.1) in Eq. (HTl) can be evaluated as 


(1.1) = az^a&z^&P™t+6(^o+6)[a;“a;^,P^Jz^a)Pi,(t'b)] = 

= aZ/a&r'bP“t^(z^a+&)P^„(z^a)K,P™,(l^&)] +ai^aKP^++,(z^a+t,)P™,(z^&)[a:a -Pmjl^a)] + 

+aVabVbP^Xh > P?n^ i^a)] [x^ , P^, (I'b)] = 

= -ai^abl^bP^tl,{l^a+b)P^^{l^a)P^^-l{l^b - 1 ) ( 1 - 1 . 1 ) - aVahVbP^JlS'^a+b)P^^-x{Va - P)PtS^b) ( 1 - 1 - 2 ) + 

+avMbP^+l^{Va+b)P?n^-,{l^a - 1)PI,-M - 1) (1-1-3) = 

= -aVaP?,tX>,i'^<^+b)PmS^a)bmbPtSvb) (1 -1 -1-H*) - aZ^aP™++, (z^a+t, )P^ )P^,-1 (z^b) (1-1-2.H*)- 

-^^Z.^™!^(^a+6)^^,(^Z.)«^aP“„(z/,)(l.l.l.P*) - hVbP!;;tlS^-+b)PlS'^b)P^..M (1-1-2-P*) + 
+am„6m,P“t',(z^a+6)P“„(z.a)P^^(z.,) (1.1.3.H*) + amaP“+^^(z.,+,)P“ Ji/a)P^,_i(z.,) (1.1.3.P*) + 

+&^b^"mt^K+6)^m„-i(^a)P^,(z/6) (l-l-3.C*) + P“+^,(z.,+,)P“_i(i/,)P^^_i(i/b) (1-1-3.P*). (42) 

The next stage of evaluation will be performed with respect to the sums over a and b after the first stage calculations 
of other terms will be completed. 

The second term can be evaluated as 


(1.2) = avMbP^i.S'^a - l)P™,(z^f, - l)[a;++b,P“+^,(z/a+b)] = 

“ Prrta+b-'^^^°-+b) 0 ‘{xna + l)6(m{, + l)Pm„ + l (z^o)Pm6 + l (z^b) (1.2.1*) + 

+P:n 1 Xb-l^^a+b)a{ma + 1)P“ „ + l (z/a)P^, (z^b) (1-2.2*) + 
+ P“t+.-l(^a+b)&(mb + l)P?nS^a)Pl, + M) (1.2.3*) + (z/„+b)P“„ (z/a)P^, (z^b) (1-2.4*). 


(43) 


The contribution (1.2.1*) describes the desirable action of the operator product j/^ vt+b te™ on the set m (cf. Eq. 

dSl)- 

Special consideration should be given to the case a = b. The involved commutator can be expressed as 

(2) = - i)[x~^xta, P^Sl'a)P^^Jty2a)] = 

= a^l2a{>^a - l)P^l(z/2a)[x-^P“ Jz/„)](2.1) +a2b.„(z., - 1)P“ Jz/„ - 2) [x+, P^^^ (z^^a)] (2.2). (44) 

The first term (2.1) in Eq. (l4^ can be expressed as 

(2.1) = - l)PmKz^2a)[a:■^P“ Jz^a)] = 

= P^yMa)a^Va{ya - 1) (-2P“ ^ _i (z., - 1) + P^^_^{l^a - 2)) = 

= -2p2“ Jb/2a)a2(z/a - l)m,P“ Jz.,) (2.1.1*) - 
-2p2“Jb.2a)az.,P“_i(:/a)(2.1.2) + 2p2“Jj/2a)aP“_i(z.,) (2.1.3*) + 
PP^yS^2a)a^ma{ma - l)P^^{Va) (2.1.4*) + P^“J:y2a)2a(TOa - l)P^„_i(z^a) (2.1.5*) + 

+P™l(^2a)P“„_2(z^a) (2.1.6*). 


(45) 
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The term (2.1.2) needs further evaluation that can be performed using identities Eq. (l40l) as 

(2.1.2) = -2a^,p2ajz.2a)P“„_l(^^a) = 

-2amaP“ Ji/a)P^l(i^2a) (2.1.2.1*) - 
-2{ama + 1 - a)Pf^^_MPXS^ 2 a) (2.1.2.2*) - 

-2P^^_MP^.S^2a) (2.1.2.3*). (46) 

The second term can be evaluated as 

(2.2) = - 1)P“ - 2)[4a,^ml(^2a)] = 

= {a^irua + 2){ma + l)Pm„+ 2 ('^a) + 2a{ma + l)Pm„+l(j^o) + 

+i"™J^a))P^l_l(^2a) = 

a^inia + 2){ma + l)Pm„+2(j^a)Pm“a-l(^2a) (2.2.1*) + 

+2a(m, + l)P“„+i(z2,)P^“_i(z.2a) (2.2.2*) + 

+P“Jz.a))p2a_i(z.2a) (2.2.3*); (47) 


Two more types of commutators (3) and (4) related to the second part of the Hamiltonian need to be evaluated. 
For the first type (a b) one has 

(3) = {a + b)i2a+b[x+x+x-^f„P^Jiya)P!^,MP^tl,{i^a+b)] = 

= {a + &)i/a+bP^+^^(!2a+&)[a;+a;+,P“ Jz/a)P^,(i^b)] (3.1) + 

+ {a + b)i2a+bP^^{i^a)P^^{i^b)[x~^b^P^+l^{iya+b)] (3.2) + 

+(a + &)i/a+b[a;+a:+,P“j!/a)P^,(!2f,)][a;“+f„P^+^^(!/a+&)] (3.3). (48) 

Each part can be evaluated separately. Consider them following the order in Eq. (|48)) 

(3.1) = {a + b)l2a+bPii,tl^{l^a+b)[xtx'l,P^^{l2a)P^^it2b)] = 

= {a + b)l'a+bPm^_^_^{l^a+b)Pma,-li’^a)Pmb-l(^b) (3.1.1) + 

+(« + b)i^,+bP^+liiya+b)Pfi.Jiya)PL,-iM (3.1.2*) + 

+ {a + b):2,+bP:ktli'^a+b)P:k.-MP^,i>yb) (3.1.3*). (49) 

Using Eq. dini) one can switch to the desirable set of basis functions for the problematic term (3.1.1) as 

(3.1.1) = {a + b)l'a+bPmt+b(’^°'+b)Pma-l(b'a)Pmi,-l(^b) = 

= Pma,-l(^a)Pmb-l(b'b){o, + b){ma+b + l)Pm('(^t + l (^^a+b) (3.1.1.1*) + 

+^m„-i(^a)Pi,_i(i^b)((a + b)ma+b + 1)P“+^,(b^a+b) (3.1.1.2*) + 

+i"™„-i(i^a)P^,_i(*2b)PX-i(^-+'>) (3-1.1-3*). (50) 


The next contribution (3.2) in Eq. (l48l) can be evaluated as 


(3.2) = (a + 6)l/a+bP^„(j^a)P^,(b'b)K+b,P^+^,(l^a+b)] = 
= -(a + &)l/a+bP^„(l^a)Pi,(l^b)P^+^,_l(l^a+b- 1) = - (o + &)ma+bP^„ (j^a)Pi, (j^b)P^++, (j^a+b) (3.2.1*) - 

-^™„(^a)P^,(i^b)P“t^,_l(*2.+b) (3.2.2*). 


( 51 ) 
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Similarly one can evaluate the remaining contribution (3.3) 

(3.3) = {a + b)va+b[xtx'^,P^Jiya)PL,M][x~^b^Pfi,+l^iiya+b)] = 

— b)^a+bPm^^^-l{>^a+b — 1) (-Pm^ (^a)^m6-l + ^ma-l (^“)Pm6 (^b) + Pm„-1 (^o)Pmb-l (^b)) = 

= -(a + &)ma+6P^+^,(^a+b)Pm„(i^a)P™,-l(t'fc) (3.3.1.1*) - 

(3.3.1.2*) - 

-{a + &)ma+bP“+^,(:/a+b)P“_i(^^a)Pi,(^^b) (3.3.2.1*) - 
-P™:!,-l(^a+b)P“„_l(^a)Pi,(^b) (3.3.2.2*) - 
-(a + &)ma+bP^+^^(z/a+&)P^„_l(z/a)P™,-l(l^b)(3.3.3.1*) - 

-P“t^,_l(j^a+b)P“„_i(j^a)P^,_l(i^b)(3.3.3.2*). (52) 

The remaining commutator (4) is similar to Eq. (l48l) but describes the case a = b. The related equation takes the 
form 

(4) = 2a:/2a[a;+^a;^a.PmJi^a)P™Kj^2a)] = 

= 2az/2aP™l(j^2a)[a:+a;+,P“_^(z/a)](4.1) + 

+ (2a):.2aP“„ (^2a)](4.2) + 

+ (2a)r/2a[x+a:+,P“Jr/a)][a:2-,P^l(^2a)](4.3). (53) 

Each part can be evaluated separately. Consider them following the order in Eq. (l53)) 

(4.1) = 2a:/2aP™“„(j^2a)[a;+X+,P“ Jz/a)] = 

= 2aZ.2aP^l(i^2a)P“„_2(^a) (4.1.1) + 

+2 ■ 2az.2aP^l(j^2a)P“„_i(j^a) (4.1.2*). (54) 

Using Eq. (l40l) one can switch to the desirable set of basis functions for the term (4.1.1) in Eq. (1541) as 

(4.1.1) = 2aZ.2aP^l(i^2a)P“„_2(i^a) = 

= P“„_2(^a)2a(m2a + l)P^l + l(^^2a) (4.1.1.1*) + 

+ P™„-2(^a)(2am2a + l)P^l(^^2a) (4.1.1.2*) + 

+Pm.-2(^a)P^l_l(*^2a) (4.1.1.3*). (55) 

The next contribution (4.2) in Eq. (l53l) can be evaluated as 

(4.2) = 2aZ.2aP“J^^a)[x2- ,P^l(z^2a)] = 

= -2aZ.2aP“„(*^a)P^l_l(*^2a " 1) = -2am2aP“ J^a)P^l(^^2a) (4.2.1*) - 

-P“Jz.a)p2“_i(z.2a) (4.2.2*). (56) 

Similarly one can evaluate the remaining contribution (4.3) 

(4.3) = 2aZ^2a[a;+X+,P“^(z^a)]fe>PmKi^2a)] = 

= -2aZ/2aP^l_l(^^2a - 1) (2P“ (z/„) + = 

= -2am2aP^l(i^2a)P“„_2(i^a) (4.3.1.1*) - 
-P^l_l(i^2a)P“„_2(i^a) (4.3.1.2*) - 
-4aTO2aP™“„(i"2a)P^„_l(j^a) (4.3.2.1*) - 

-2p2“_l(z.2a)P“„_l(^^a) (4.3.2.2*). (57) 


Finally one have to express all results as sums. Unchanged products of Laguerre polynomials are skipped in each 
sum for the sake of simplicity. We begin with Eq. (1431) and introduce new notations given after each result in its final 
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form designated with the letter F. 


^ + + + + {Fl.l.A)' 

b,a^b 

(1.2.2 + ^ + l)P“„+i(^a) = 

b,a^b 

= Y P™.-i(^c)a(m, + l)P“^+i(;/J (Fl.l.P)-^p2a_i(i/2a)a(m, + l)P“^+i(i/J (Fl.l.P.l); 


c,a<c a 

(1.2.4*)sum = 2 ^ -frnl+b-l('^a+f>) = 

b.a^b 

= i^(a-l)P“_i(z.,) (i^l.l.C)-i^p2a_i(^.2j (Fl.l.C.l). 

a a 


(58) 


The next group of expressions Eq. (HU) can be evaluated as 


(1.1.1.A + B*)sum = - Y, ^^abtUb = -nm + Y, o^^ama {F1.2.A)] 

b,a^b a 

{1A.2.A + B*)sum = - Y = 

a,b^a 

= -n^P“_i(:/a) (P1.2.P.1) + 

a 

+ ^az.,P“_i(z.J (1.2.P.2); 

a 

(1.2.P.2) = ^az/aP“„_i(i^a) = (P1.2.P.2.A) + (P1.2.P.2.P) + (P1.2.P.2.C'); 

a 

{F1.2.B.2.A) = ^oma = m {F1.2.B.2.A); 

a 

(P1.2.P.2.P) = ^(ama + 1 - a)P“„_i(i^a) (P1.2.P.2.P); 

a 

(P1.2.P.2.C) = (P1.2.P.2.C); 

a 

{lAA.A*)sum = ^Y (^'bnabmb = ^ (P1.2.C'.l) - ^ X! (PI.2.(7.2); 

b.a^b a 

(1.1.3.P + (7*)s„m = Y “”^“-f’mi,-l(^h) = 

b,a^b 

= mYP?n^-M (P1.3.7l)-^amaP“„_i(z^a) (P1.3.P); 

a a 

(lA.'i.D*)sum = 2 ^Y, Pma-lib'a)Pmb-lib'b) (Pl-4). 


(59) 
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Next we consider the contribution of the case a = b from Eq. (1441) for the first part of the Hamiltonian 


(2.2.1*)sum — 2 ® + 2)(TOa + l)f’m„+ 2 (^a)^m 2 a-l(^ 2 a) (^2.1); 

a 

(2.2.2*)sum = Cl{ma + l)^’m„ + l(i^a)^’m 2 a-l(^ 2 a) {F2.2)] 

a 

{2.2.3*)sum = 2 (^a))^m 2 n,-l(^ 2 a) (^"2.3); 

a 

(2.1.1*)™™ = -J2a^(i,^- l)ma (F2.4); 

a 

(2.1.2.1*)™™ = - awg = -m (F2.5.1); 

a 

(2.1.2.2*)™™ = - y^(awa + 1 - a)P™_^_i(r/a) (F2.5.2); 

O 

(2.1.2.3*)™™ = - (F2.5.3); 

O 

(2.1.3*)™™ = (P2.6); 

a 

(2.1.4*)™™ = ^ X! ^"^^airria - 1) (P2.7); 

a 

(2.1.5*)™™ = i ^ 2a(ma - l)P“„_i(z^a) (P2.8); 

a 

(2.1.6*)™™= i^p2“Jr^2a)P“„_2K) (i"2.9); 

a 


( 60 ) 
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The next contribution is associated with the remaining part of the Hamiltonian (a 7 ^ b) Eq. (I48|) 

= 2 + b){ma+b + ^)Pm^^^ + i{l^a+b) (^ 3 . 1 ); 

b,a^b 

(3.1.1.2*)s„j„ = 2 Pma-li'^a)Pm,i-l{b'b){{Oi + b)ma+b + ^) (E3.2.1); 

b,a^b 

(3.1.1.3*)s„m = 2 ^ma-l(^a)^m 6 -l(^ 6 )^ml+i,-l(*^a+&) (E3.2.2); 

(3.1.2 + 3*)su-m = X! + b)l/a+bPmh-l(^b) = 

a,b^a 

= (F3.2.3.H) - ^ 2ai.2,P“_i(z.,) (F3.2.3.P) - 

a a 

(3.2.3.C')- ^ (3.2.3.P); 

o a,c<.a 

(3.2.3.C) = (P3.2.3.C'.l) + (P3.2.3.C'.2) + (P3.2.3.C'.3); 
(P3.2.3.C'.l) = - ^oTOa = -m (P3.2.3.C'.l); 

a 

(P3.2.3.C.2) = - ^(ama + 1 - a)P“^_i(j^a) (P3.2.3.C.2); 

a 

(P3.2.3.C'.3) = -J2p^^_^{ua) (P3.2.3.C.3); 

a 

(3.2.3.P) = - ^ (:.,) = 

a,c<a 

- J2 c(m, + l)P;(,^+i(i.,)P“ (P3.2.3.P.1)- 

o,c<a 

- ^ (cm,+ l)P;(,^(i.,)P“ (P3.2.3.P.2)- 

a,c<a 

- ^ P^,-l(^^c)P“„_l(^^a) (P3.2.3.P.3); 

a,c<a 

(3.2.1*)s„m = -i ^ (a + b)ma+b = “^ X! “"^o(a - 1) (P3.2.4) + ^ X! (P3.2.4.H); 

b,a^b a a 

{3.2.2*) sum = — 2 = 

b,a^b 

= (^3.2.5)+ i^p2“_i(z.2,) (P3.2.5.H); 

a a 

(3.3.1.1*) 

sum + (3.3.2.1*) 

sum — E aWaP^,_i(i^&) (P3.2.6) + 

b,a{>b 

+ y^2am2aP^^_i(t^a) (P3.2.6.H); 

a 

(3.3.1.2*)sum + (3.3.2.2*)sum = ~ = 

b,a^b 

= -\Y. Pm.-MP^s-M) (P3.2.7) + 
+ E^™l-l(^ 2 a)P“_l(^a) (P3.2.7.H); 

a 

(3.3.3.!*)™™ =-i ^ (a + 6 )ma+bP“^_i(z.a)P^^_i(i/,) (P3.2.8); 

b,a^b 

(3.3.3.2*)s„™ = —— P™)'(^^_i(i^a+&)Pm^_l(t'a)P™^_l(^'b) (P3.2.9). 

b,a^b 


( 61 ) 
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Finally we consider the second part of the Hamiltonian for a 
as 


b. The contributions from this part can be expressed 


+ l)P™l+i(z^2a) (F4.1); 

a 

(4.1.1.2*),„„ = + (P4.2.H); 

a 

(4.1.1.3*),,„= i^P“_2K)P^l_i(z^2a) (F4.2.B); 

a 

{4.1.2*Um = (F4.3). 

a 

{4:.2.1*)sum = y^2aTO2a (P4.4.H); 

a 

(4.2.2*U^ = - i^p2a_^(z.2,) (F4.4.B); 

a 

(4.3.1.1*)sum = -i^2am2aPma-2('^a) (P4.5.H); 

a 

{4.3.1.2*)sum = -iE^-l-i(^ 2 a)P“„_ 2 K) (F4.5.B); 

a 

{A.3.2.U),um = -^^4am2aP“„_i(i^a) (F4.5.C); 

a 

{A.3.2.2*)sum = (F4.5.P). 


(62) 


Finally we collect all contributions with like terms. We begin with the terms corresponding to the same or mod- 
ihed sequences {m} representing the partitions of the same number m (e. g. (F4.1) or (F4.4.A)), while the terms 
representing the partitions of different numbers (e. g. (F4.5.A)) should cancel each other. Indeed, the contribution 
(F4.1) is represented by the sequences with modified parameters TOq —>■ — 2, m 2 a —>■ W 2 a + 1 conserving the sum 
m = aViia, the contribution (F4.4.A) is represented by the same sequence (diagonal term), while the contribution 
(F4.5.A) is represented by the sequences modified as rua —>■ rna — 2 representing the partition of the different number 
m — 2a. 

First we collect all diagonal terms entering with the same basis function V'{m}({'^})- These terms are enumerated 
by indices (F1.2.A), (F1.2.B.2.A), (F1.2.C.1), (F1.2.C.2), (F2.4), (F2.5.1), (F2.7), (F3.2.3.C.1), (F3.2.4), (F3.2.4.A), 
(F4.4.A). The sum of all terms can be evaluated as 


2 ^ 

— nm + arua (F1.2.A) +m (F1.2.P.2.A) + ^ (^1.2.67.1) - a^ml (P1.2.C.2) - 

a a 

- ^ a^ma{va - 1) (i^2.4) - m (F2.5.1) + ^ E “ 1) (-^2.7) - m (F3.2.3.C'.l) - 

a a 

^E amaia - 1) (F3.2.4) + ^ E 2 «”^2a (P3.2.4.A) - E 2 a»^2a (P4.4.A) = -nm + — 11. 


(63) 


This expression taken together with the first term n(n— l)/2 in Eq. (l34l) coincide with the diagonal term in Eq. 

Next we collect the off-diagonal terms represented by the partitions of the same number m. They include contri- 
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butions (Fl.l.A), (F2.1), (F3.1) and (F4.1) which can be written as 


= + + (Fl.l.A); 


b,a^b 


nf. = + + (F2.1); 

a 

lf± = IY1 ^“„-i(^a)P^,_i(^^;,)(a + 5)(m,+, + l)F“+^^+i(z.,+,) (F3.1); 

ir-- = ^E^-.-2K)2a(™2a + l)P^l + l(^^2a) (F4.1). 


b,a^b 


(64) 


The notations reflect the related anharmonic interaction actions for processes (a, b) i - > (a+6) in separately considered 

cases (a ^ b) 01 a = b. 

The contributions Eqs. (l63l) and (l64l) added together are equivalent to the recursive equation Eq. (1321) . It is obvious 
for the diagonal term. Eor off-diagonal terms it becomes clear if we consider the total expression for the wavefunction 
amplitude Eq. (1291) which will appear in the right hand side of Eq. p2p with several polynomial products modified 
as in Eq. ([M]). For instance the contribution (Fl.l.A) can be expressed in the form 


- E 


E ‘^{rnyytVb Va+b ({:/})] V'/ 

\{m} 




(65) 


Since the action of the operator product VaybUl+b expression in parenthesis Eq. (1651) does not change this 

expression because it shifts the summation indices without modifying the whole sum (the “border” terms with rua = 0 
or rrih = 0 do not contribute because of the factor matrif,) one can rewrite this sum shifting the summation indices in 
the internal sum as ruo —>■ rria + 1, mj, —>■ mb + 1, ma+b —t 'ma+b — 1 in the form 


E*w(M)5E rriambya, % 2/j+6d{m}V’/ 

. {m} b,a^b 


( 66 ) 


The expression multiplied by the basis function 'h{m}({^}) is identical to the first term in the right hand side of 
Eq. (15^ for different indices a and b. Similar transformations can be applied to the other three terms in Eq. (IMl) 
reproducing all other contributions in Eq. (15^ which can be then obtained setting the coefficients with the like terms 
'h{m}({j^}) equal to zero, provided that all other off-diagonal terms corresponding to the partitions of a numbers 
different from m are compensated with each other. One should notice that although the solution for the amplitude 
djmj. satisfying Eq. (1321) if nontrivial is definitely the eigenstate of the problem other solutions can exist because the 
basis functions 4'|^}({z/}) are not necessarily linearly independent. For instance for n = 3, m = 2 the basis functions 
tp{ 2 ,o) = “ l )/2 — + 1/2 and V'(o,i) = Pi {^ 2 ) = 1^2 — 1/2 are related to each other as ^/’( 2 ,o) = ~'0(o,i) 

for all integer partitions of n = 3 including vz) = (3, 0,0), (1,1,0), (0, 0,1). This nonorthogonality is a source 

of trivial solutions. However, according to numerical study all eigenstates probed yet for n < 25 can be represented 
by integer number sequences as described in the main body of the manuscript. 

Consider the remaining off-diagonal terms, which indeed cancel each other as one can see below. First we can perform 
some obvious cancellations including (F1.2.F.2.F) -|- (F2.5.2) = 0, (F1.2.F.1) -|- (F3.2.3.A) = 0, (F3.2.2) -|- (F3.2.9) = 
0, (FI.2.5.2.17) -H (F3.2.3.17.3) = 0, (Fl.l.C) -h (F3.2.5) = 0, (Fl.l.C.l) -h (F3.2.5.A) = 0, (FI.1.5.1) -h (F2.2) = 0, 
(F1.4) -h (F3.2.1) -h (F3.2.3.5.3) -h (F3.2.7) -h (F3.2.8) = 0, (F4.2.5) -h (F4.5.5) = 0, (F3.2.6.A) -h (F4.5.17) = 0, 
(F4.3H(F3.2.3.5) = 0, (F1.1.5)-H(F3.2.3.5.1) = 0, (F2.3)-f (F4.4.5) = 0, (F2.5.3)-h(F2.9)-f (F4.2.AH(F4.5.A) = 
0, (FI.3.5) -I- (F2.6) -I- (F2.8) = 0, (F3.2.7.A) -I- (F4.5.5) = 0. The sum of remaining terms is also equal to zero as 


shown below 
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(Fl.S.A) - J2{ama + 1 - a)P^^_M (F3.2.3.C.2) - 

a a 

- Y, {afna + l)P:i,MP^.-M {F3.2.3.D.2)- ^ cm,P^^_,{ua) (F3.2.6) = 


c,a<c 


c,a<c 


a 

+E^-.- 


, - ^ cmc - E 


CTTlc — arUa 


c:a<c 


c:a>c 


+ 


—arua + arria + E 1 — (a — 1) 


= 0 . 


(67) 


Thus the validity of Eq. (15^ and, correspondingly, of Eq. dMl) is proved. 


GENERATING EUNCTION FOR WAVEFUNCTION AMPLITUDES 

Since eigenfunction amplitudes for a certain principal quantum number n can be expressed through the eigenfunc¬ 
tions of the same problem with different principal quantum number m Eq. (I29|) one can continue this procedure 
adding more numbers until the last quantum number approaches zero. Then possible eigenfunctions of the problem 
with the principal quantum number n can be described by a sequence of numbers = n,ki,k 2 , ...kp-i^kp = 0. 
The wavefunction for a certain integer partition {v} of number n representing the resonant sequence of popu¬ 
lation numbers satisfying the conservation law Eq. (l30l) can be expanded over the basis of the products of Laguerre 
polynomials i^{m}iW}) as 


^ ^ ) ( 68 ) 
{m} 

where the sequences {m} represent all integer partitions of the number ki. This procedure can be written in the 
matrix form 


c = M"’'=id, (69) 

where “vectors” c and d represent wavefunction amplitudes in the partition spaces for numbers n and fci (m) and the 
matrix element between partitions of numbers n and ki is given by the corresponding product of Laguerre polynomials 
Eq. (l28l) . Remember that if fci = 0 one can set all amplitudes dy^y = 1. 

Repeating the procedure Eq. (1551) until reaching the last non-zero number in the sequence {A:} one can express the 
eigenfunctions in terms of the product of matrices M by the vector of all unities Ifcp_i of the size equal to the number 
of integer partitions of the number kp-i as 


■p-2 

. 2=1 


(70) 


This expression is complicated and we cannot evaluate it in general. Yet it is possible to calculate the related 
generating function defined as 


G{.{v},{y})= E Y yp 

{m}i {”i }2 


V{m}2({"i}i)- E yp-i 

{m}p_i 


(71) 


where r]({m}) = is a number whose partition is realized by a sequence {m}. The wavefunction amplitudes 

for the specific partition {i/} and sequence {k} can be found as the coefficients of the generating function Eq. (1711) 
polynomial expansion with the products 

The generating function can be calculated using the identity for Laguerre polynomials which reads (29| 


m—O 


( 72 ) 
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Using this identity one can evaluate the generating function Eq. dzlD. First this function can be separated into the 
products of contributions of Laguerre polynomials of different arguments 1/j where j stands for the phonon state. 
Then the most right sum in Eq. m for the specific j can be evaluated using Eq. as 


Sp-iU) = exp - 


yi-i 


{l + y: 






(73) 


where m{p — 2, j) is the value of number in the sequence {m}p- 2 - The next to the most right summation in Eq. 
(1711) yields 


c / o I yp-2 yp—2yp—i yp—i\ ^ i j i j j a,?) 

Sp- 2 {j) = exp I —--^ I (1 + j /^_2 + y;_ 22 /p_i) 


(74) 


Repeating this procedure p — 1 times we can express the contribution of the specific state j as 


S{j) =exp I -- 


E n 

/.ff(l</<S<p) *=/ 


VI 


1 + 2/i + {yiy2y + ...{yi..yp-iy 


(75) 


where is the population number of the photon in the state (partition) {ly} of interest. 

Finally taking the product of contributions Sj over all states j = 1,2.. Eq. (ESI) and using the identity !3 = 

— ln(l — x) we obtain the generating function in the final form 


G(M,{y})= n 

/.s(l</<s<p) 


1 - 




n 

n[i 


■yf 


(yiy2y + ...{yi-yp- 


(76) 


used in the main body of the manuscript. 

Eq. (17^ is too complicated for general analysis of the problem eigenstates. Yet it is sufficient to fully characterize 
the eigenstates containing non-zero contribution of the single phonon state \n > represented by the population number 
sequence = 5kn- We restrict the consideration to strictly decreasing sequences {fc} just because we know from 
the solution that all eigenstates of interest can be found using these sequences. Then for the state Vk = 5kn one can 
leave only unity from the second product in the right hand side of Eq. (1761) because otherwise the power of some of 
the variables yi will approach or exceed n which is not acceptable since these powers should form strictly decreasing 
sequence n > ki > k 2 > ks... > kp = 0. The only acceptable choice of the contribution from the first product satisfying 
the sequence selection requirements can be made taking the composite product of yi x (2/12/2) x (2/12/22/3)x (2/i2/2---2/p-i)- 
The related sequence is given by {n,p—l,p—2,p—3, ...0) and it defines the eigenstate 4'(p). There are n such sequences 
determined by the integer number p changing from 1 to n. In all cases the absolute value of the amplitude of the 
wavefunction in the single phonon state is equal unity. 

To illustrate the results of the previous paragraph consider the case of p = 3. Then the part of the generating 
function that is of interest is given by 


g 

P{M,yi,y2)= n ^ 

/.s(l</<s<p) _ *=/ 

= (1 - 2/i)(l - 2/2)(l - 2/12/2) = 1 - 2/1 - 2/2 - 2/12/2 + 2/12/2 + 2/i2/2 + 2/12/2 “ 2/?2/2- (77) 

It is clear that the only contribution of interest is associated with the term 2/12/2 since the power of 2/1 should be larger 
than the power of 2/2 and they should be both strictly positive. This is in a full accord with the previous conclusion. 

One can prove using the generating function methods (cf. Eq. (IM1) 1 that the wavefunctions 'I'(p) are normalized 
by 1 (see Sec. ,). Using the normalization factor as in Eq. p3p we find the normalized by 1 wavefunction amplitude 
for the state |n > to be c„ = lj\/n. Since we found n states containing the single phonon state and n|c„p = 1 this 
gives an additional evidence than our choice is complete. This result is used within the main body of the manuscript 
to investigate the time evolution of the probability for the system to be in a single phonon state \n >. 
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NORMALIZATION OF THE WAVEFUNCTION «'(p) 

The proof is separated into two parts. In the first part we prove that the longest possible strictly decreasing sequence 
(n,n — 1, ...0) results in the wave function 4'(n) with amplitudes defined as 

n 

CM = (-1)^'= (-1)”, n = ^ k,.,. (78) 

k=l 

The normalization of this wavefunction by 1 is proved in Sec. □ because the eigenstate considered there differs from 
Eq. (1781) by signs only. In the second part the proof of the normalization by one will be given for the function (p) 
defined by the general sequence {n,p — l,p — 2, ...0) {p < n). 

To prove the first statement assume that it is valid for p = n — 1 and consider its generalization for p = n. Then, 
according to Eq. (1681) we can express the wavefunction amplitude for the partitions {v} of the number n in terms of 
the wavefunction amplitudes for the partitions {m} of the number n — 1 as 

/ 

cM = (-ir-En (79) 

{m} k 

The sum is taken only over the sequences {m} representing the partitions of n — 1. This expression can be 
evaluated using the generating function M{{v},x) defined as 

M{{iy},x) = (-l)”-^En (80) 

{m} k 

where the sum is taken over all sequences {m} of nonnegative population numbers. The wavefunction amplitude of 
interest is given by the expansion coefficient of the function M{{i/},x) with the factor 

One can express the generating function Eq. (I79p as a product of independent sums for all integer fc’s and then 
evaluate each sum using Eq. (|72l). Then we get 

MiM, x) = tiZE (1 _ . (81) 

k 

Assume that I is the minimum index k in the product in Eq. (1811) corresponding to the non-zero population number 
i//. Then using the algebraic identity (1 — x’‘)/{l — x) = X]g=o can represent Eq. (|5T|) as 

1 — 1 n 

M(M,x) = (-1)"-iE*M1 -^T'”' n (82) 

q=0 k=l+l 

The term of interest with the power is given by the highest power of x term in Eq. (1821) which reads 

(_l)"a;"-i(_l)Efc^'., (83) 

The coefficient with this term has the form of Eq. (1751) . which proves the first statement. 

Consider the proof of the second statement. According to Eqs. (1551) and (1751) the amplitudes of the wavefunction 
'k(p) can be expressed as 


CM = (-l)”-'E^^-^'”''^(lA)(-l)'”^ (84) 

{m} k 


where the summation (^^) is taking over integer partitions of a number p (^j. kruk = p).The normalization of the 
corresponding wavefunction is given by (cf. Eq. (1231) 1 


E 


\cm\ 


n 






M ^^*=1 *■ {m} \ k / {r} \ k / 


(85) 
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where the summations is over partitions of number n kvk = n) and number p kruk = p, = P)- 

As previously the summations over partitions in Eq. (1851) can be evaluated using the generating function method 
introducing the new generating function 

R{x,y^,y2) = I E ( HE . (86) 

{m}\ k / {r}\ k ) 


The normalization factors of interest Eq. (1851) can be found as the expansion terms accompanying the products 

n P P 

x"yiy2- 

The generating function can be calculated similarly to Eq. (EU evaluating the sums over specific states k as 

y^ + yl + xHi-yi){i-y2y 


Sk = exp 


(87) 


and then taking the product of all these expressions for k = 1,2.... Then the generating function can be expressed as 

(1 - xyi){l - xy 2 ) 


R{.x, j/i, j/2) ■ 

It is convenient to reexpress this generating function as 

R{x,yi,y2) + (1 _ 2,)(i _ xyiy2)' 


( 88 ) 


(89) 


The first term does not contribute to the terms of interest (x^'y^y^ with n > p) while for the second all expansion 
coefficients for n > 1 are equal unity, which proves the normalization of the states 'I'(p) by 1. 


SOLUTIONS FOR PRINCIPAL QUANTUM NUMBERS n < 5 

Below we describe the construction of the basis of eigenstates of the anharmonic Hamiltonian V Eq. (1221) for 
principal quantum numbers n < 5 to illustrate the algorithm proposed within the main body of the manuscript. We 
use the strictly decreasing sequences of quantum numbers {A:} satisfying the additional constraint ki-i—kt > ki — ki+i. 
One can easily check that all obtained solutions correspond to the wavefunctions normalized by 1 and orthogonal to 
each other and all obtained energies are consistent with Eq. 6 in the main text defining energies as a function of 
generating sequence in the form 


e{{k}) = 


n(ji — 1) 
2 


p-i 

+E “ 1) “ kih+i]. 


(90) 


i=0 

In addition to constructing the eigenstate basis we also discuss the symmetries of the states and their overlap with a 
single phonon state derived in the main text. 


n=l 

Here the basis consists of only one partition of the number 1 which is = 1. The anharmonic Hamiltonian is 
represented by a zero diagonal element. The only available sequence in this case is (1,0) corresponding to all unity 
solutions = 1 and zero energy in agreement with Eq. (TOl) . 


n=2 


Here the basis consists of two partitions {vi = 0, 1/2 
within this basis as 


1) and (2,0). The operator h Eq. (l27t can be represented 


h2 = 



(91) 
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There are two generating sequences including (2,0) and (2,1, 0). The first one is represented by all unity eigenvector 

, while the second one should be expressed using the Laguerre polynomials basis set with the second quantum 

number m = 1. Since there is only one partition for m = 1 (see Sec. ^ the solution can be expressed in terms of 

the single Laguerre polynomial ^^(1) = — 1 = ^ , while the coefficient with this polynomial has to be one 

because of the final step 1 —>■ 0 corresponding to the all ones solution. The energies of these states are given by ±1, 
respectively, in agreement with the analytical prediction of Eq. ([90]). 

Two eigenstates can be expressed through each other using the transformation e —>■ —e, Cjj,} —>■ Cjj,} • (— 
described in the main text. Both states contain the single phonon state (0,1) with the amplitude absolute value equal 
to 1 in agreement with the derivation in the main text. 


n =3 


Here the basis consists of three partitions (0, 0,1), (1,1, 0) and (3,0, 0). The operator h Eq. ([27|l can be represented 
within this basis as 

/O 3 0\ 

^3 = 2 0 1 (92) 

\0 3 0/ 


There are three generating sequences satisfying the algorithm proposed within the main text including (3,0), (3,1, 0) 
and (3,2,1,0). The first one generates the stationary state given by the all ones eigenvector 1 , and the second 




one generates the state given by (1) = z/i — 1 = 0 similarly to Sec. □ The third one should be expressed 

2 

/ 1/2 \ 

by the superposition of Laguerre polynomials ^^(1) = — l)/2 — + 1/2 = I —1/2 1 (partition (2,0) for 

VI/ 2 ; 

/-l/2\ 

Pi = 2) and ^^(1/2) = 1^2 — 1/2 = ( 1/2 1 (partition (0,1) for pi = 2). The coefficients with these polynomials 

V- 1 / 2 ; 

are determined by the remaining subsequence (3,2,1,0)_ = (2,1,0) corresponding to the stationary state ^ 

(coefficients vi — 1 for the partitions of pi = 2, see Sec. 0 ) leading to the final expression for the wavefunction 
amplitudes 


^^(3.2,1.0) = - L^-'^(l/2) = I -1 1 ■ 


(93) 


The energies of the associated stationary states are given by £( 3 , 0 ) = 3, £( 3 , 1 , 0 ) = 0; £( 3 , 2 , 1 , 0 ) = ~3 in agreement with 
the analytical prediction of Eq. dSI). 

Pair of eigenstates described by the sequences (3, 0) and (3, 2,1, 0) can be expressed through each other using the 
transformation £ —>■ —£, —>■ ■ (— 1 )^*''' described in the main text while the state determined by the sequence 

(3,1,0) transfers to itself. All three states contain the single phonon state (0, 0,1) with the amplitude absolute value 
equal to 1 in agreement with the derivation in the main text. 
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n=4 


Here the basis consists of five partitions including (0,0,0,1), (0, 2,0, 0), (1,0,1,0), (2,1,0, 0) and (4, 0,0, 0). The 
operator h Eq. (EH) can be represented within this basis as 


h4 


/O 2 4 0 0\ 

4 0 0 2 0 

3 0 0 3 0 

0 14 0 1 

\0 0 0 6 0 / 


(94) 


There are five generating sequences including (4,0), (4,1,0), (4, 2, 0), (4, 2,1,0) and (4, 3, 2,1, 0). The first sequence 
corresponds to all ones eigenstate with energy e 4 _o = 6. The second sequence corresponds to the state 


l/’(4.1,0) — L 


(I'l-b 


(1) = z/1 - 1 = 


/- 1 \ 

-1 

0 

1 

V3y 


(95) 


with energy e 4 _i^o = 2. The third sequence creates the state corresponding to the symmetric combination of two 
Laguerre polynomials because the coefficient with polynomials do suppose to be equal one for the terminating step 
2—^0. This state can be expressed as 


'0(4,2,O) — ^2 ^(1)+-^!^ ^(1/2) — — l)/2 — 1^1 + 1^2 — 


2 

-1 

0 


V2; 


(96) 


with energy €420 =0. The fourth sequence corresponds to anti-symmetric combination of the same functions as in 
Eq. (IMl) formed similarly to Eq. (IMl) as 


V’(4,2,l,0) — -^ 2 "^^ ^^(1) — ^^(1/2) — — l)/2 — — 7^2 + 1 


-1 

0 

-1 




(97) 


with energy 64 , 2 , 1,0 = ~2. The eigenstate determined by the remaining sequence (4, 3, 2,1,0) can be expressed using 
the algorithm Eq. (|29)) and the reduced sequence solution ■i/'( 3 , 2 ,i,o) Eq. (|M1) as 


V'(4.3,2,1.0) — '0(3,2,1,O)(1)E3 ^ (1) -|-^(3,2,l,0) (2)Ti ^ '*(1)L)^ ''(1/2) -|-■!/7(3, 2 , 1 , 0 ) (3)Ti ^ "'(1/3) — 

-1 ^ ^ ^ -1)(.. - 2) j _ ^ ^ 


-(!^1-1)/,n r(l'2-l)/ 


-(‘^3-1) 


(98) 


The energy corresponding to this sequence is given by —6. 

All energies agree with the theory predictions Eq. dmi). Pairs of eigenstates described by sequences (4, 0) and 
(4, 3, 2,1,0) and (4,1,0), (4, 2,1,0) can be expressed through each other using the transformation e —>• —e, —>■ 

C{j,} • (—l)^**"’ described in the main text while the state determined by the sequence (4,2,0) transfers to itself. 
Four states determined by the sequences (4,0), (4,1,0), (4, 2,1,0) and (4, 3, 2,1,0) contain the single phonon state 
(0, 0,0,1) with the amplitude absolute values equal to 1 in agreement with the derivation in the main text. 
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n =5 


Here the basis consists of seven partitions (0,0,0,0,1), (0,1,1,0,0), (1,0,0,1,0), (1,2,0,0,0), 
(3,1,0, 0,0) and (5,0, 0,0, 0). The operator h Eq. (E7)) can be represented within this basis as 


( 2 , 0 , 1 , 0 , 0 ), 


/O 5 5 

6 0 0 
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4 0 0 
0 4 4 
0 1 6 
0 0 0 


0 0 
3 1 

2 4 
0 0 
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3 6 


o\ 

0 

0 

0 

0 
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(99) 


\0 0 0 0 0 10 0 / 


The sequences corresponding to the eigenstate basis are given by (5,0), (5,1,0), (5,2,0), (5,2,1,0), (5,3,1,0), 
(5, 3, 2,1,0) and (5,4, 3, 2,1,0). Corresponding eigenstates calculated using the algorithm Eq. (|29ll similarly to 
the previously considered cases can be expressed as 
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The corresponding energies of stationary states are given by 10, 5, 2, 0, 

-2,- 

-5, 

-10. 


All energies agree with the theory predictions Eq. (ITOl) . Pairs of eigenstates described by sequences (5,0) and 
(5,4,3, 2,1, 0), (5,1,0) and (5, 3, 2,1,0), and (5,2,0) and (5, 3,1,0) can be expressed through each other using the 
transformation e —>■ —e, • (—1)^« described in the main text while the state determined by the sequence 

(5, 2,1,0) transfers to itself. Five states determined by the sequences (5,0), (5,1,0), (5, 2,1,0), (5, 3, 2,1,0) and 
(5,4,3, 2,1, 0) contain the single phonon state (0, 0,0,1) with the amplitude absolute value equal to 1 in agreement 
with the derivation in the main text. 


MATLAB PROGRAMS TO CALCULATE HAMILTONIAN, EIGENSTATES AND EIGENENERGIES OF 

THE FPU a PROBLEM 

Below we describe the supplied Matlab functions targeted to calculate eigenstates and eigenenergies of the problem 
Eq. (ITfl) using the proposed algorithm Eq. dZOl). The verification of the theory using these functions is straightforward. 
For instance consider the quantum number sequence (5, 3, 1, 0). The eigenfunction and eigenenergy corresponding 
to this sequence can be found numerically using the command line call “[v,E,vn] = EigSt{[5, 3, 1, 0]);”. The 
outcomes include the eigenvector v of the modified problem Eq. (EH), the normalized by 1 eigenvector vn and the 
eigenenergy E = —2 all calculated using the recursive algorithm based on Eq. (|7ni). The result can be verified 
generating the system Hamiltonian as “y = InitAnhHInfN{5)-’’. Then the standard Hamiltonian “H” can be found 
using “H = y.E[". The operation ‘‘‘‘H*vn-E*vn" should return a vector of zeros (with the appropriate accuracy which 
is 10“^® in my computer; the equation is typed in the way that it can be copied and pasted directly to the Matlab 
command line at least in my computer). The modified representation of the Hamiltonian can be found using the 
command “yl = InitAnhElInf NMod{b)\”. Then the operation “yl.H*v-E*v” should lead to a vector of zeros as 
well. The reader is strongly encouraged to try these codes for different quantum number sequences to examine the 
validity of the proposed solution. 

Before running the commands you need to copy all supplied Matlab files to the same folder and change the current 
folder used by Matlab to that specific folder. 


Function y = FindResModes(M, S, N) 

This function finds all integer partitions of a number S made using N — M + 1 integer numbers M, M + 1, ...N. It 
is targeted to find the basis set for the problem of interest expressed in terms of phonon population numbers. The 
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call “y = FindResModes{l, N, N)” returns all integer partitions of the number N. Partitions are expressed as rows 
of the matrix representing the answer. 

The number S must be less or equal N and the parameter M must be greater or equal to 1 and less or equal to 
both N and S. This function should work fine until N < 60. 

This function is independent of other author’s programs. 


Function y = InitAnhHInf N{N) 

This function generates resonant Hamiltonian of the third order anharmonic interactions within the long-wavelength 
and large size limits for the principal quantum number N. 

The outcome y.H returns the Hamiltonian matrix, the other outcome y.Hsp returns the same matrix in the sparse 
matrix Matlab form (type “y = InitAnhHInfN{5); y.H” in the Matlab command line to generate the Hamiltonian 
for N = 5). 

The number N must be an integer number greater than 1 and less or equal to 30; for TV > 30 the memory can be 
insufficient for the resulting Hamiltonian matrix. Then the only sparse matrix outcome can be used. 

This function depends on the function ’’FindResModes” generating the basis set of partitions. 


Function y = InitAnhHInf NMod{N) 

This function generates modified Hamiltonian h (Eq. (1271) 1 of the third order anharmonic interactions within the 
long-wavelength and large size limits for the principal quantum number N. 

The outcomes y.H returns the Hamiltonian matrix, the other outcome y.Hsp returns the same matrix in the sparse 
matrix Matlab form (type “y = InitAnhHInf NMod{5); y.H” in the Matlab command line to generate the modified 
Hamiltonian for = 5) 

The number N must be an integer number greater than 1 and less or equal to 30; for TV > 30 the memory can be 
insufficient for the resulting Hamiltonian matrix. Then the only sparse matrix outcome can be used. 

This function depends on the function ’’FindResModes” generating the basis set of partitions. 


Function [y, E, yn] = EigSt{B) 

This function generates the eigenstate and the related eigenenergy using the sequence of quantum numbers B 
determining this state. 

The vector y returns the eigenstate amplitudes in the modified population number representation and the number 
E returns the eigenstate energy, while yn returns the wavefunction amplitudes normalized by 1 in the true population 
number representation. 

The set B must be a strictly decreasing set of integer numbers beginning with the principal quantum number n 
and ending with 0 (For example eigenstate and eigenenergy corresponding to the set (5,3, 2,1,0) can be found typing 
^\y,E] = EigSt[[f>, 3, 2, 1, 0])”). 

This function uses the functions “FindResModes” to generate partitions, “BasFun” to make the wavefunction 
expansion over the basis of the polynomial products and ’’NormFact” to switch to the standard basis. 


Function y = CollectsEigsPartit{N) 

This function is expected to generate the full set of eigenstates and the related eigenenergies using the special set of 
sequences of strictly decreasing quantum numbers with the constraint that the array of differences for any sequences 
is non-decreasing (the set (5,3, 2,1, 0) is acceptable, while the set (5,3, 0) is not because 5 — 3 < 3 — 0). 

The matrix y.V returns eigenvectors in a modified population number representation as columns, the row y.Etst 
returns eigenenergies corresponding to eigenvectors and the matrix y.Conib returns the matrix of sequences used to 
generate the eigenstates. 

The completeness of the basis set can be tested calculating the rank of the matrix of eigenvectors “A = rank{y.V);”. 
The result can then be compared with the size of the matrix of eigenvectors that can be determined using the command 
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“S' = rankijj.V);" . According to our consideration the rank of the matrix is identical to both sizes of that matrix at 
least up to = 25. Unfortunately, we cannot prove the completeness of the basis set analytically for arbitrarily N. 
N must be a positive integer number. 

This function uses the functions “FindResModes” to generate partitions and “EigSt” to find eigenstates and energies 
for each sequence. 


Function y = NormFact{P) 

This function generates the normalization factors for conversion between actual (normalized by 1) and modified 
basis sets. 

The vector y returns the column of normalization factors for each row of the input matrix P. 

Matrix P is a matrix of sets (rows) of population numbers representing some basis states. All numbers should be 
integer and nonnegative. 

This function is independent of other author’s programs. 


Function y = BasFun{Part, Bas) 


This function calculates the basis functions as a products of Laguerre polynomials Eq. (1281) . The sequence Part 
represents the sequence of population numbers (partition) and the matrix Bas represents the sets of orders of corre¬ 
sponding Laguerre polynomials which can be expressed as the set of partitions. 

Vector y returns the vector of the products of polynomials of the specific population number set Part for each 

row of the sequence Bas. For instance if the function arguments are the matrix Bas = f ^ ^ partition 


Part = (5, 2,1) then the outcome vector will he y = f 

All sequences should contain only integer nonnegative numbers. 

This function uses the function ’’Coeft” evaluating Laguerre polynomials. 


0 0 1 


Function y = Coeft{r, k, n) 

This function calculates the basis function for the individual state k with the population number n expressed in 
terms of associated Laguerre polynomials [1^ as L^~'^\l/k). 

Inputs r and n must be nonnegative integer numbers. 



